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Abstract 



The gamma distribution with small shape parameter is difficult to characterize. 
In particular, standard algorithms for sampling from such a distribution often fail, 
^ ' so something special is needed. In this paper, we first obtain a limiting distribution 

for a suitably normalized gamma distribution when the shape parameter tends to 
00 I zero. Then this limiting distribution provides insight to the construction of a new, 

simple, and highly efficient acceptance-rejection algorithm. Pseudo-code, and an 
CN . R implementation, for this new sampling algorithm are provided. 

^ . Keywords and phrases: Acceptance-rejection method; asymptotic distribution; 

characteristic function; gamma function; R software. 

^ ■ 1 Introduction 

.9.' 

Let y be a positive gamma distributed random variable with shape parameter a > 0, 
denoted by F ~ Gam(a, 1). The probability density function for Y is given by 

= y>0, 
T{a) 

where the normalizing constant, T{a) = y'^~^e~y dy, is the gamma function evalu- 
ated at a. This is an important distribution in statistics and probability modeling. In 
fact, since the gamma distribution is closely tied to so many important distributions, 
including normal, Poisson, exponential, chi-square, F, beta, and Dirichlet, one could ar- 
gue that it is one of the most fundamental. But despite its importance, when the shape 
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parameter a is small, the gamma distribution is not so well understood. Consequently, 
seemingly routine calculations can become intractable. For example, in the popular R 
software, the functions related to the gamma distribution — in particular, the rgamma 
function for sampling — become re latively inaccurate when the shape parameter is small 



fiR Development Core Teamll201ll . "GammaDist" documentation). To circumvent these 
difficulties, and to move towards new and more accurate and efficient software, an im- 
portant question is if Gam(a, 1), suitable normalized, has a meaningful non-degenerate 
limiting distribution as a — )■ 0. In this paper, we give an affirmative answer to this ques- 
tion, and then we use the result to develop a new and improved algorithm — in terms of 
both accuracy and efficiency — for sampling from a small-shape gamma distribution. 

When the shape parameter a is large, it follows from the infinite-divi sibility of the 
gamma distribution and Lindeberg's central limit theorem ( lBillingsleylll995l . Sec. 27) that 
the distribution of Y is approximately normal. Specifically, as a — )■ oo. 
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a, 



— N(0, 1) in distribution. 



For large finite a, better normal approximations can be obtained by working on different 
scales, such as log Y or Y^^^. Our interest is in the opposite extreme case, where the shape 
parameter a is small, approaching zero. Here, neither infinite-divisibility nor the central 
limit theorem provide any help. In Theorem [T] below, we prove that —a\ogY converges 
in distribution to Exp(l), the unit-rate exponential distribution, as a — t- 0. Our proof 
is based on a basic property of the gamma function and Levy's continuity theorem for 
characteristic functions. 

Motivated by the limit distribution result in Theorem [H we turn to the problem of sim- 
ulating from a small-sha pe gamrna distr ibuti on. This i s a cha llenging problem with many 
proposed solutions; see Devrove ( 1986 ) and Tanizakj ( 2008 ). For sr nall shape parame- 
ters, t he d efault methods implemented in R and MATLAB, due to lAhrens and Dieter 
(119741 ) and iMarsaglia and TsangI (l2000h , respectively, have some shortcomings in terms 
of accuracy and/or efficiency. The exponential limit in Theorem [T] for the normalized 
gamma distribution suggests a convenient and tight envelope function to be used in an 
acceptance-rejection sampler (IFlurvi Il990l ). We fiesh out the details of this new algo- 
rithm in Section [3] and provide an R function to implement this method online. This new 
method is simple and very efficient. Compared to the method implemented in rgamma, 
which frequently returns logF = — oo when a < 0.001, ours can quickly easily produce 
genuine samples of log Y for even smal ler a. It i s also more efficient than the ratio-of- 
uniform samplers recently proposed in IXi et al.l (120131 ). Our proposed strategy should 
also be useful for other problems. 



2 Limit distribution result 

We start with some notation. For the gamma function T{z) defined above, write f{z) = 
logr(z). Then the digamma and trigamma functions are defined as fi{z) = f'{z) and 
f2{z) = f"{z), the first and second derivatives of the log gamma function f{z). Recall 
that these are related to the mean and variance of logy, with Y ~ Gam(a, 1): 

E,(logF)=/i(a) and V„(logF) = /2(a). 
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These formulae are most directly seen by appl ying those we ll-known formulae for means 
and variances in regular exponential families f lBrownlll986l . Corollary 2.3). Next, write 
Z = —a log Y. To get some intuition for why multiplication by a is the right normaliza- 
tion, consider the foUowiiig recu rrence relations for the digamma and trigamma functions 
(lAbramowitz and Stegunlll966l . Chap. 6): 

/i(a) = /i(« + !)-!/« and /2(a) = /2(« + 1) + 

Then, as a -> 0, 

E,(Z) = -afi{a) = -ah{a + 1) + 1 = 0(1), 
V„(Z) = a'hia) = a^Ma + 1) + 1 = 0(1). 

That is, multiplication by a stabilizes the first and second moments of logy. Towards a 
formal look at the limiting distribution of Z, define the characteristic function 



(1) 



ip^it) = E«(e**^) = E^y-^"*) = T{a - tat)/T{a), 
where i = is the complex unit. 

Theorem 1. For Y ~ Gam(a, 1), —alogY — Exp(l) in distribution as a — j- 0. 

Proof. Set Z = — alogF. The gamma function satisfies T{z) = T{z + l)/z, so the 
characteristic function fait) for Z in ([T]) can be re-expressed as 



(fait) 



T{a — iat) T{1 + a — iat)/{a — iat) 



1 Til 



r(«) 



r(l -t-a)/a 



1-it r(l + o„) 



where Oa are terms that vanish as a — )■ 0. Since the gamma function is continuous at 1, 
the limit of (fa{t) as « — exists and is given by 1/(1 — it). This limit is exactly the 
characteristic function of Exp(l), so the claim follows by Levy's continuity theorem. □ 



3 Simulating small-shape gamma variates 



Simulating random variables is an important problem in statistics and other quantitative 
sciences. For example, simulation studies are often used to compare performance of com- 
peting statistical methods, and the Bayesian approach to statistical inference relies heavily 
on simulations since the posterior distributions on moderate- to high-dimensional param- 
eter spaces rarely have a tractable closed-form expression. Here we shall demonstrate 
that the limiting distribution result in Theorem [1] helps provide an improved algorithm 
fo r simulating gamma rand om variables with small shape parameter, compared to that 
of Ahrens and Dieterl ( 1974 ) implemented in the rgamma function in R. 

For Y ~ Gam(a, 1) with a near zero, let Z = —alogZ. To simulate from the 
distrib u tion o f Z, one might c onsider an acceptanc e -rejec tion scheme; see, for example. 



Langd fll999l . Chap. 20.4) or iGivens and Hoetind (120051 . Chap. 6.2.3). For this, one 



needs an envelop function that bounds the target density and, when properly normalized, 
corresponds to the density function of a distribution that is easy to simulate from. By 
Theorem [1] we know that Z is approximately Exp(l) for a ^ 0. More precisely, the 
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density ha{z) of Z has a shape hke for z > 0. Therefore, we expect that a function 
proportional to an Exp(l) density will provide a tight upper bound on ha{z) for z > 0. 
We shall similarly try to bound ha{z) by an oppositely-oriented exponential-type density 
for z < 0, as is standard in such problems. 

The particular bounding envelop function rjaiz) is chosen to be as tight an upper 
bound as possible. This is done by picking optimal points of tangency with ha{z). For 
this, we shall need a formula for ha{z), up to norming constant, which is easily found: 

ha{z) = ce~^~^ ''"^ , z e (—00,00). 
The normin g constant c satisfies = r(a + 1). By following standard techniques, as 



described in iLangd (119991 . Chap. 20.4), we obtain the optimal envelope function 



ce ^, for 2 > 0, 
cw\e^^ , for 2; < 0, 



where A = A(a) = — 1 and w = w{a) = a/e{l — a). A plot of the (un-normalized) 
target density ha{z) along with the optimal envelope ria{z) is shown in the left panel of 
Figure [U for a = 0.1. The normalized envelope function ria{z) corresponds to the density 
function of a mixture of two (oppositely-oriented) exponential distributions, i.e., 

1 -Exp(l) + -^{-Exp(A)}, 



1 + w 1 + w 

which is easy to sample from using standard tools, such as runif in R. 

Pseudo-code for a new program, named rgamma.ss, for simulating a sample of size 
n from a small-shape gamma distribution, on the log scale, based on this acceptance- 
rejection scheme is presented in Algorithm [1] R code for this program is also available 



at the first author's website (www.math.uic.edu/~rgmartin). Though the R code is 



already fast, a proper implementa tion should be done in C, with inline functions, test 
sample recycling (e.g. RubinI 1976 ). and perhaps other optimizations. 



The acceptance rate r{a) for the proposed method is 

r{a) = {1 + w{a)}-' = {l + -^^^}"\ (2) 

which is plotted in the right panel of Figure [1] as a function of a near zero. It is clear 
from the graph, and also from the approximation r{a) ~ 1 — a/e for a ^ 0, that the 
acceptance rate converges to 1 as a — )■ 0. This indicates the high efficiency of the proposed 
acceptance-rejection method when a is small. This is to be expected based on Theorem[Tl 
when a ^ 0, Z is approximately Exp(l), so an algorithm that proposes Exp(l) samples 
with probability 1/(1 + w) pa 1 will almost always accept the proposal. 

One can easily experiment with the new function rgamma. ss and the function rgamma 
built in to R. As the R documentation for rgamma indicates, for a < 0.001, one will 
frequently observe Y = 0, corresponding to samples logF = —00. So, obviously, that 
function is not suitable for small a values. On the other hand, the new function rgamma . ss 
can easily produce genuine samples of log Y for even smaller a values. Moreover, the code 
is very simple and the sampling is highly efficient. 
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Algorithm 1 - Pseudo-code for program rgamma. ss to simulate n independent samples 
of logy, with Y ~ Gam(Q;, 1) and small shape parameter a. 
set A -h- A(a), w ^ w{a), and r ^ r(a) as in the text, 
for z = 1, . . . , n do 
loop 

f/^Unif ^Unif denotes the random number generator 

if t/ < r then 

z i \og{U /r) 

else 

z ^ log(Unif)/A 
end if 

if ha{z) / ii]a{z) > Unif then 
Zi-i^ z 
break 
end if 
end loop 

log^i i Zi/a 

end for 
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Figure 1: Left: the (un-normalized) ha{z) and envelope ria{z) for a = 0.1; Right: the 
acceptance rate r{a) in ([2]) as a function of a. 
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